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We have investigated the six most tightly bound states of each of the atoms, helium and lithium, in 
intense magnetic fields. In this regard, we report new data for two new states of lithium that have not 
been studied thus far in the hterature. The energy levels of these first few low-lying states are calcu- 
lated in this study employing a pseudospectral method as the computational tool. The methodology 
involves computing the eigenvalues and eigenvectors of the generalized two-dimensional Hartree-Fock 
partial differential equations for these two- and three-electron systems in a self-consistent manner. 
The method described herein is applicable to calculations of atomic structure in magnetic fields of 
strength in the so-called intense-field regime, as it exploits the natural symmetries of the problem 
without assumptions of any basis functions for expressing the wave functions of the electrons or 
the commonly employed adiabatic approximation. It is seen that the results obtained herein are 
improvements upon single-configuration calculations as well as configuration-interaction methods. 
It is also seen that the pseudospectral method employed here is considerably more economical, from 
a computational point of view, than previously employed methods such as a finite-element based 
approach. 
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I. INTRODUCTION 

The motivation to study atoms in magnetic fields of 
strength beyond the perturbative regime is attributable 
in a large part, due to the discovery of such fields being 
present in white dwarf stars [UUj and neutron stars [HIS]. 
The most commonly observed neutron stars - pulsars, 
have been observed to have magnetic fields on the order 
of 10^ - lO^T [5]. Magnetars |7], which are strongly mag- 
netized neutron stars, can have magnetic field strengths 
well in excess of lO^T. White dwarf stars on the other 
hand have somewhat less extreme fields, albeit still high, 
~ 10^ - lO'^T [B]. At such high field strengths, a Zeeman- 
type perturbative treatment of the field \E is not possible. 
The structure of atoms is considerably altered from the 
low field case. 

The problem of atoms in magnetic fields has been tack- 
led by various researchers since the 1970's using a variety 
of different methods. In the literature, there exist numer- 
ous studies of hydrogen [9Hl8] and many recent studies 
of helium [TM57] atoms in strong magnetic fields. There 
have also been studies conducted for molecules and chains 
of atoms for both hydrogen and helium atoms in strong 
to intense magnetic fields [38U45J . Moreover, our recent 
investigation |46] using single-configuration Hartree-Fock 
(HF) theory [57] was seen to yield accurate upper bounds 
for the binding energies of hydrogen and helium in strong 
magnetic fields. Our later study [^, obtained accu- 
rate binding energies for helium and lithium atoms in 
strong magnetic fields. The numerical method employed 
in the latter study was seen to yield considerable divi- 
dends, in terms of reducing the computational time for 
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the solution of the problem. The atomic structure code 
provided therein is readily portable and could easily be 
included in existing atmosphere models of neutron stars 
and white dwarf stars, without much computational ex- 
pense. Presently, the reader is referred to Ref. [4^ for a 
short chronological review of both one- and two-electron 
systems in strong magnetic fields. 

In sharp contrast to these somewhat simpler systems, 
there is very limited work available in the literature for 
atoms with more than two electrons. A brief chrono- 
logical review of the topic is provided herein. One of 
the first studies to investigate atoms in intense magnetic 
fields, in particular the iron atom, was by Flowers et al 
[3S] in 1977. This variational study extended the work 
due to the authors in Ref. ^5Dj and obtained binding en- 
ergies of iron atoms and condensed matter in magnetic 
fields relevant to neutron stars. Errors in this study were 
later corrected by Muller [5T]. Other methods included 
density functional studies [EH [S3] and also employed the 
Thomas-Fermi-Dirac method [54j|55] for estimating bind- 
ing energies of atoms in intense magnetic fields. The first 
comprehensive HF studies of atoms with more than two 
electrons were carried out by Neuhauser et al (SU E7j for 
magnetic fields greater than lO^T, thus being directly rel- 
evant to neutron stars. Elsewhere, HF studies of atoms 
and molecules in intense magnetic fields were conducted 
by Demuer et al |58j . with results consistent with previ- 
ous findings. All of the above treatises, Refs. [49H58] . con- 
cern themselves with magnetic fields in excess of lO^T, 
well into the so-called intense magnetic field regime. At 
these field strengths, the interaction of the electron with 
the nucleus of the atom becomes progressively less dom- 
inant, in comparison to its interaction with the field it- 
self. One of the first studies to carry out a rigorous 
HF treatment of atoms with more than two electrons 
in strong or intermediate field strengths was Ref. [22]. 
Therein, they obtained estimates of the binding energies 



of a few low-lying states of lithium and carbon atoms, 
in low to strong magnetic fields. Elsewhere, Ivanov |59j 
and Ivanov & Schmelcher [551 ISHl IS0HS3] , have over re- 
cent years, carried out detailed HF and post-HF studies 
of atoms with more than two electrons using a numerical 
mesh-method for solving the unrestricted HF equations 
[55] , The special meshes were so constructed as to fa- 
cilitate finite-difference calculations in a two-dimensional 
domain using carefully selected mesh node points [M] , 
They were able to ascertain the binding energies of the 
first few low-lying states of low-Z atoms such as lithium, 
beryllium and mid-Z atoms such as boron and carbon 
etc., using this method. Moreover, using a gaussian basis 
of functions for expressing the wave functions of the elec- 
trons [2SH31], adopting a full configuration-interaction 
method, Al-Hujaj & Schmelcher [65j .66, have been able 
to estimate the binding energies of lithium and beryllium 
atoms in strong or intermediate magnetic fields, thereby 
improving upon previously obtained results. The sodium 
atom in a strong magnetic field has also been studied by 
Gonzalez-Ferez & Schmelcher [67] with estimates of bind- 
ing energies for the same. Elsewhere, low lying states of 
the lithium atom have also been studied in strong mag- 
netic fields using a configuration-interaction method, em- 
ploying the so-called freezing full-core method both with 
[55] and without [51] correlation between electrons. Re- 
cently, Medin & Lai [HJ |32] have also studied atoms and 
molecules and infinite chains of condensed matter in mag- 
netic fields greater than lO^T, using density-functional- 
theory. Mori et al [351 [IZ] have studied mid-Z atoms 
in strong to intense magnetic fields using perturbation 
theory as well, obtaining results consistent with previ- 
ous findings. Finally, as mentioned earlier, recent work 
on the lithium atom using a single-configuration HF cal- 
culation employing spectral methods for solution yielded 
accurate eigenvalues and eigenvectors of low-lying states 
in strong or intermediate magnetic field strengths [48] . 
The method was seen to be computationally more eco- 
nomical than conventional HF solution methods using 



either finite-element techniques 
methods [6]. 



46J or finite difference 



In the literature, there are very few studies of atoms 
with more than two electrons in strong magnetic fields. 
The methods developed for the determination of bind- 
ing energies are either very involved such as Quan- 
tum Monte-Carlo and configuration-interaction tech- 
niques using different basis functions, or they are com- 
putationally demanding such as numerical mesh-methods 
and finite-element techniques. Configuration-interaction 
(CI) methods rely upon determining the expansion coeffi- 
cients and assume that the individual electron wave func- 
tions are already pre-determined from essentially several 
single-configuration calculations, i.e., the orbitals are not 
optimized during the calculation. A variant of post-HF 
methods that optimizes the orbitals in addition to the ex- 
pansion coefficients is the so-called Multi-Configuration 
Hartree-Fock (MCHF) method 70J. Both these post-HF 
methods yield considerable improvements with regard to 



the estimates of the upper bounds for the energies of 
various states. In the intermediate range of magnetic 
field strengths, where both the nucleus of the atom and 
the magnetic field have interactions with the electrons 
that are approximately equal in magnitude, the single 
configuration approximation then becomes increasingly 
ineffective with greater number of electrons. However, 
these methods are computationally more intensive than 
a single configuration calculation. Thus far, the most ac- 
curate CI methods involve decomposing the wave func- 
tions into a Gaussian basis set relying upon separation 
of variables in cylindrical coordinates [e.g. |551|55]. These 
methods do however require a large set of basis func- 
tions. On the other hand, MCHF methods would require 
fewer basis functions, as the orbitals get optimized during 
the computation with the coefficients [e.g. iTOj. The sep- 
aration of variables and/or basis decompositions speed 
up the computation in these post-HF methods consid- 
erably. However, there do not exist hitherto, any fully 
two-dimensional (2D) post-HF studies of multi-electron 
atoms in intense magnetic fields. This is partly due to 
the computational overhead associated with adopting a 
fully two-dimensional picture. Central to the develop- 
ment of such a method would be the fast and accurate 
computation of the single-configuration problem in a full 
2D framework without any basis expansions and separa- 
tion of variables. Wave functions so determined could be 
used directly in 2D configuration-interaction calculations 
or the problem could be cast into a MCHF framework. 
In the latter, the orbitals within each configuration could 
be optimized as the calculation proceeds. Such post-HF 
studies would yield more accurate data for the structure 
of atoms in intense magnetic fields. The aim of the cur- 
rent study is to facilitate this larger goal, by providing a 
fast and accurate solution of the full 2D single configu- 
ration problem of atoms in intense magnetic fields. 

The method outlined in the current study is an exten- 
sion of the method developed in our two previous studies 
Refs. [351 HH]. The framework does not make any as- 
sumptions of basis functions and neither is it restricted 
to the adiabatic approximation. The computations are 
seen to be readily extendable to arbitrary field strengths 
and atoms with greater number of electrons. The over- 
all algorithm of solution using pseudospectral methods is 
also computationally straightforward to implement and 
has already been seen to yield considerable dividends in 
terms of computational time [48j . Obtaining accurate es- 
timates of the energy levels of atoms, in particular low-Z 
atoms, in strong and intense magnetic fields will ulti- 
mately facilitate a proper understanding of the spectra 
of neutron stars and white dwarf stars. Thus the central 
aim of the current work is to provide accurate estimates 
of the binding energies of the first few low-lying states of 
the simplest low-Z atoms; helium and lithium, in strong 
and intense magnetic fields using an unrestricted two 
dimensional single-configuration pseudospectral method. 
Our subsidiary goal is to additionally provide a computa- 
tionally efficient and economical software that could be 



directly integrated into atmosphere models of neutron 
stars. 



II. THE HF EQUATIONS 

We shall begin with the generalised single- 
configuration HF equations for an atom with 'rig' 
electrons and nuclear charge 'Z', in a magnetic field 



that is oriented along the z-direction. A derivation of 
the single-configuration HF equations can be found in 
our earlier work [35], here we shall present only the 
salient points. The single configuration HF equation 
can be written in cylindrical coordinates as shown in 
Eq. [T] Here the length scale is in units of Bohr radii 
and the energy is scaled in units of Rydberg energy 
in the Coulomb potential of charge Ze (see below for 
definitions). 
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where i, j = 1, 2, 3, ..., rig and r^ = \J pf + zf . The total Hartree-Fock energy of the state is given by 
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The direct and exchange interactions are determined according to the method outlined in Ref. [46], as the solutions 
of the elliptic partial differential equations for the potentials given by. 



V2$c = -4^|V,(p„z,)|' 
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and 
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aE{pi,Zi) = -Aii%p*Api, Zi)il)i{pi, Zi). 
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The wave function of a given configuration of electronic 
orbitals is assumed to be given by a single Slater deter- 
minant as. 



$ = A„ 



[V'l,V'2,V'3,---,'0ne-l,V'r, 



(5) 



where A„^ is the anti-symmetrization operator. The in- 
dividual electronic wave functions V'i are given by. 



Vi. = ^,U,z.)e''"*'X.(5»), 



(6) 



where i labels each of the rie electrons. The two- 
dimensional single particle wave functions ipi{pi,Zi) are 
taken to be real functions. 

The integration with respect to the azimuthal coordi- 
nate, (j), has been carried out, prior to writing the result 
in Eq. ([T]) above. The contribution due to electron spin 
has also been averaged out a priori. It is to be men- 
tioned in this regard, that in the current study we shall 
only be concerned with fully spin-polarised states (FSP), 
in other words all the electrons of the atom are assumed 
to be anti-aligned with the magnetic field. Such states 
have an exchange interaction between the electrons pro- 
viding an extra coupling term in the HF equations, as- 
Additionally, the FSP states are the seen to be the most 
bound states in increasing magnetic fields in the intense 
field regime. The extension to partially spin-polarised 



configurations is easily achieved by eliminating the ex- 
change term in the HF equations. In the current study, 
we have chosen to work in units of Bohr radii along with 
the definitions given below. 

The Bohr radius for an atom of nuclear charge Z is 
given by as/Z , where as — h/ arrive is the Bohr radius of 
the hydrogen atom. The magnetic field strength parame- 
ter Pz^ is given by the expression [3z — B/{Z^Bo), where 
Bq is the critical field strength at which point the transi- 
tion to the intense magnetic field regime occurs [6J . This 
is defined as Bq = {2a'^m^'^c^)/{eh) w 4.70108 x lO^T. 
Thus, beyond a value oi (iz ~ 1, the transition to the 
intense magnetic field regime occurs and the interaction 
of the electron with the nucleus becomes progressively 
less dominant as Pz increases. The energy parameter 
of the j*'' electron is defined as e^ = Ei/{Z'^Eao), with 
EoD = hct^rn^c^, the Rydberg energy of the hydrogen 
atom. For brevity we shall refer to the units of energy as 
Ez,oo, which should be remembered as the Rydberg en- 
ergy in the Coulomb potential of charge Ze. The quantity 
a — e^/(47reo^c) w 1/137 is the fine structure constant. 
In the current study, all the physical constants were used 
in SI units. Additionally, the magnetic field B, is taken 
to be in units of Tesla. The above written Eq. (fTl) rep- 
resents the A^-coupled Hartree-Fock equations in partial 
differential form for an A^-electron system with nuclear 



charge Z. The system of equations is solved iteratively; 
see Section ITVl for the numerical details. 

Based upon the above definition of f^z-, it is convenient 
to classify the field strength [7T] as low {jiz < 10"'^), 
intermediate, also called strong (10~^ < /3z < 1) and 
intense or high (1 < /^^ < oo). It is immediately evident 
upon inspection that Eq. (Il]) is a set of coupled non- 
linear partial differential equations. The equations are 
coupled through the exchange interaction term between 
the electrons. Presently, we shall describe in detail, the 
numerical methodology developed in the current study. 



III. 



THE PSEUDOSPECTRAL APPROACH 



The HF equations collectively given by Eq. ([T]) above, 
describe a coupled eigenvalue problem. The atomic state 
function consists of a single Slater determinant that de- 
scribes the wave function. The numerical solution of the 
coupled equations in Eq. ([I]) proceeds via the so-called 
self consistent field (SCF) method due to Hartree [47] . 

First we find a solution to the hydrogenic problem; 
Eq. (fTl) without the direct and exchange interactions. 
This allows us to obtain an initial estimate for the wave 
functions. Second, using these estimates the elliptic par- 
tial differential equations for the direct and exchange in- 
teraction potentials in Eqs. ([3| and Q are solved. With 
these potentials now obtained, the coupled HF problem 
including the direct and exchange interactions in Eq. ([l| , 
is solved as a full non-linear eigensystem. The eigen- 
values obtained are the individual particle energies e.; 
and the normalized eigenvectors are the wave functions, 
4'iiPijZi)- The SCF iterations then proceed with the up- 
dated electron wave functions and the steps from the sec- 
ond step described above, are repeated until convergence. 

Central to the entire scheme is the method of solution 
of the partial differential equations. For a numerical so- 
lution of the above scheme, we employ a discretization 
based on a pseudospectral approach as described in an 
earlier paper [48] . In contrast to our earlier work, in the 
current study we employ a cylindrical coordinate system. 
As a result the methodology for setting up the problem 



is considerably different from that described in Ref. |48) . 
This section is arranged as follows. First we describe 
the methodology employed for solving the hydrogen atom 
using pseudospectral methods in cylindrical coordinates. 
Particular emphasis is placed on the implementation of 
boundary conditions. Thereafter, the treatment is ex- 
tended to the particular case of the helium atom in a 
single configuration and a generalization of the scheme is 
then provided for greater than two electrons. 
A. The Hydrogenic Problem 

We begin with the Hamiltonian for the hydrogenic 
problem (single-electron) in a strong magnetic field, 
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i'l {pi,Zi) ^ e^^Xi {pi, Zi) . (7) 



The solution of the eigenvalue problem in Eq. (l7| yields 
the individual electron energies and wave functions in a 
given configuration. The wave functions thus obtained 
form the initial guess to begin the Hartree-Fock itera- 
tions. The domains of both the radial and the axial co- 
ordinates are Q ^ p,z < oo. The problem maintains 
azimuthal symmetry and thus, a solution of Eq. (l7| in 
this domain, when reflected about the z = plane (re- 
specting z— parity of course) and revolved about the z- 
axis through 27r, gives the solution in three-dimensional 
cylindrical coordinates. 

First the semi-infinite domains are compactificd using 
the transformations, 



and 



2 tanh (p^) — 1 



i/i = 2 tanh(zi) - 1. 



(8) 



(9) 



This ensures that the semi-infinite domain is mapped 
from [0,(X))(8)[0, cxd) to [-1, 1]® [-1, 1]. With the transfor- 
mations given in Eqs. (Sl) and Q, we can re- write Eq. (It]) 
as 
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-1 f^ + x 



^{3-y^-2y)il + y)^ + im,Y 



-^-Al--A^-y'-^y)% 



tanh 



1 -l-a; 



dx 4 

-I -2 



tanh 



1 + 2: 



[tanh-^(i±^)]V[tanh-i(lf^)]' 



A (x, y) = CiVi (a;, y) ,(10) 



where, we have dropped the subscripts on the coordinate labels. The discretisation points are thereafter taken to 



be the commonly used Chebyshev-Lobatto points [151 [7^ 
[73] given by 



Xj — cos {nj/N) 



(11) 



where j = 0, 1, ..., A'^. As is customary, we employ monic 
polynomials of degree N as the cardinal functions to in- 
terpolate between these points and are given by [IHl [72] , 



1 ^ 



(12) 



with 



■' k=0 



N 



n(^ 

fe=0 



Xk) 



(13) 



Derivatives of these interpolating polynomials at the dis- 
cretisation points then yield the so-called Chebyshev dif- 
ferntiation matrix, whose elements are given by, 



A,= 



1 ^ 

■' fe=0 
k^i,j 



Xk) 



^j\Xi ^k) 



i^^j) (14) 



To illustrate, let us consider the simple case of ap- 
plying a zero or null-type Dirichlet condition at the end 
points of the domain (x ^ —I, j = N and x = I, j = 0). 
This implies that we are fixing the values at the end 
points to zero, i.e., u{xk'=i) = and u{xk'=N+i) = 
in Eq. ([T7| . Therefore the first and the last equation 
in Eq. ( )17p vanish identically. Thus, the first and last 
columns of the matrix D^ can be discarded as they mul- 
tiply with u{xk'=i) = and u{xk'=N+i) = 0, respec- 
tively. Similarly, the first and last rows of the matrix D^ 
can also be discarded as we are ignoring both the first 
and the last equation in Eq. (17 1 as f{xk'=i,N+i) = 0. 



Therefore the effect of imposing Dirichlet conditions at 
the end points has the effect of trimming the working 
matrix by removing appropriate rows and columns. Let 
us define this trimmed version of the Chebyshev Differen- 
tiation matrix as Dj 
labelled /,fc' = 2,... 
of the domain Xji=2, 



k' ■ The trimmed matrix has indices 
N. Thus the solution in the interior 
,,^N can be found according to 



PiXj') 



N 

E 

k'=2 



m- 



j'k 



f{xk' 



,iV, 



(18) 



where [D^) ^ is the inverse of D^ . The vectors p{x) and 
f{x) are explicitly given by. 



and 



N 

Dj] = Y^{xj 

fc=0 



xkY 



(15) 



With the differentiation matrices thus formed, it is pos- 
sible to write down a differential equation in matrix form 
as follows. Consider the differential equation, 



dx 



^u{x) ^ fix). 



(16) 



This can be written using the Chebyshev Differentiation 
matrix as, 



P = 



Pix2) 
PiXs) 



p{xn-i) 
. p{xn) . 



f{x2) 

fixs) 



f{xN-l) 
. fixN) . 



(19) 



Notice that the vectors p and f also have indices running 
from 2 to iV since the first and last rows have been re- 
moved due to the imposition of null type Dirichlet bound- 
ary conditions. On the other hand, if the Dirichlet con- 
ditions were of the fixed variety with u{xk'=i) = a and 
u{xk'=N+i) = b, then we can write Eq. (17) as 



N+l 



Y, D}k'P{xk') = fix 



J') J 



,N + 1, 



(17) 



where D^ is the square of the matrix D and pix) is the 
polynomial interpolant approximating the solution uix). 
In Eq. (17) j' and k' refer to rows and columns of the 



matrices and vectors and hence start from 1 rather than 



as in Eq. ( 11 ). Using this matrix formulation, it is then 



possible to write more complicated differential equations 
as matrix equations [481 I72| . provided the domain has 
been suitably compactified to [—1,1]. The differential 
equation in Eq. ( 16 ) has thus been converted to a dis- 



crete equation with an operator (a matrix) acting upon 
a vector. A simple solution (u) of the above equation 
can then be obtained by multiplying both sides from the 
left with the inverse operator, (Z)^)~^, once appropriate 
boundary conditions have been imposed. 



N 

Y, D}k'Pixk') = fixk') 

k'=2 



«^/,l 



^^/,^+l 



(j' 



,iV). 



(20) 



In this case, the solution in the interior of the domain 
Xj'=2,...,N can be written as, 



N 



E 

k'=2 



Pixr)=^iD^)-\,,,ifixk-) 



aDl, 



bD},N+i) 



if = 2,...,N). 



(21) 



It is then immediately evident that if a = 6 = 0, then we 
recover Eq. (18). 



Neumann-type boundary conditions are handled in a 
similiar manner by extending Eq. (21 1 using the differ- 
entiation matrix D [74j . For example, if we wish to 



solve Eq. (17) subject to the boundary conditions that 
u{xk'=N+iT= and u'{xk'=i) = 0, then we can re-write 
Eq. ( pOl by replacing 6 = and a — p{xk'=i), the latter 
a priori; 



is unknown 



N 

Yl ^fk'Pi^k') = f{xk') - D^'.iP{xk'=i) 

fc'=2 



(/ = 2,...,iV). 



(22) 



Putting E = D^^.,, with /,fc' = 2,...,N and Eq = 
[ 



D| I, D^ 1, ■■•, Djf ]^] , we can re-cast Eq. (22 1 in matrix 



form as, 



Ep = {-Eop{xi) 



(23) 



However, since u'{xk'=i) = it implies that the deriva- 
tive of the solution at the x = 1 boundary must vanish. 
We can write this condition as 



N 



N 



y^ Di,k'P{xk') = Di^ip{xi) + ^ Di^k'P{xk') = 0. 



k' = l 



k'=2 



(24) 

Substituting Bq — -Di,i and Bi == [131,2, -Di,3, .-., i?i,Ar] 
we can re-write Eq. ([24| as 

(25) 



Bopixi) + Bip = 0. 



This allows us to obtain an expression for the unknown 
quantity p{xi) in terms of p as 



Pixi) 



-B^'Bip. 



Substituting Eq. ( 26 ) in Eq. ( 23 1 we obtain 



{E-EoB^'Bi)p^{. 



(26) 



(27) 



The solution in the interior can then be obtained as, 

p^{E-EoB^^Bi)-H. (28) 

This then yields a straightforward linear algebra problem 
of solving a system of equations with is readily handled 
using standard methods. The method outlined above 
for imposing boundary conditions and formulating the 
problem, has been extended to the case of two dimensions 
and is applied to obtain a solution of the problem in 



Eq. ( 10 ) which is described below in brief. 



At this juncture it is to be mentioned that in writing 
Eq. (10 1, we have removed the co-ordinate singularity 
at X = —1, by replacing it with x — —1-1-5, where 
6 — 10~^^ in units of Bohr radii. This approximation 
produced acceptable results within error tolerances. The 
singularity at the outer boundary of the domain, given 
by a; = 1 (corresponding to p = oo), is taken care of by 
imposing a Dirichlet boundary condition since the wave 
function must vanish at infinity (see below) . Similarly, at 
the other outer boundary; y = +1, we impose a Dirichlet 
condition since here too, the wave function should vanish 
at infinity. The remaining boundaries at x,y = —1 can 
have either Dirichlet or Neumann boundary conditions, 
depending upon the behaviour of the wave function. The 
following discussion delineates the methodology for the 
2D problem. 



1. An Explicit Example - Domain Discretization 

We consider here an explicit example to illustrate the 
use of pseudospectral methods for solving an eigenvalue 
problem. In particular, we shall delve into the method- 
ology adopted here for imposing boundary conditions in 
considerable detail. The method developed here is a non- 
trivial extension of the one developed by the authors in 
Ref. [71]. 

Let us begin with a domain [—1,1] [—1,1] which is 
discretised using A^ -I- 1 points in each of the two Carte- 
sian directions; x and y, with iV = 3. This is illustrated 
in Fig. (fTl) The partial differential equation in Eq. ( 10 1 is 




( » • • n (^.^2) ■ 

(x-2,y2) {xi,y2) M 

n > > i ^^B 

(X2.y-i) (xi.ys) (xom) ■ 

= -1 Ma: = +1 m 



FIG. 1. Pictorial representation of the domain [—1, 1]® [— 1, 1] 
discretised using A'^ -I- 1 points in each direction, with N — 3. 
The outer boundaries have Dirichlet conditions imposed. 

two-dimensional therefore, following the procedure out- 
lined in Refs. |3S1[71], we can construct two-dimensional 
operators by employing Kronecker products of matrices. 
For example, a differential equation of the form. 



can be written in matrix form as, 

Dl u{x^,yj)+bl u{x^,yj) = J{x^,yj) 
z,j = l,...,iV+l. 



(29) 



(30) 



and Dy the outer 



Note that in the matrices the D 
boundaries aX x,y = +\ have already been excised due to 
Dirichlet boundary conditions and the squares of the re- 
spective differentiation matrices have been appropriately 
trimmed. As a result, the indices i^j are limited to a 
lower value of i,j = 1 rather than i,j = 0. The oper- 
ators D^ and D^ are then given by kronecker products 



with the identity matrices, 



bl = [D,, X D^] ® ly 



and 



Dt 



[Dy X Dy] , 



(31) 



(32) 



where Ix and ly are identity matrices of dimension N^ x 
Nx and Ny x Ny respectively. In our example, the number 
of node points in the x and j/-directions are equal, thus 
Nx = Ny = N . Therefore, using the above formalism we 
can write down Eq. ( 10 ) in matrix form as. 



diag(a) x l)^ + diag(6) y. 0^+ diag(c) x bl 



+diag(d) X by + diag(e) V» = -^i^i = ^ii'i, (33) 



where D„ 



D^ 



ly and Dy 



/j. (g) Dy. Since we are 



imposing Dirichlet boundary conditions at the two outer 
boundaries, this implies that we can remove the first row 
and column from each of the matrices D^ and Dy (see 
discussion above and Refs. |1H1[72]), prior to forming the 



operators in Eqs. (31) and (32). The diagonal matrices 



a, 6, c, d and e are the coefficients of the different terms in 



Eq. ( 10 ), in the order in which they are written; explicitly. 



of values at the collocation points, we can write the ma- 
trix as an extended vector comprising of the different 
columns, one followed by another. This then forms an 
7V^ X 1 vector rather than an iV x A^ matrix. Explicitly, 
we can reshape the matrix p{xi ,yj), with values given at 
the collocation points (see Fig. IT]), from 



pixi,yi) p(x2,Vi) p{x2,yi) 

P{X1,V2) P{X2,V2) p{xz,V2) 

p{.xi,yz) p{x2,Vi) p{x3,Vi.) 



(35) 



to 



'p{Xi 


2/1)1 


p{xi 


2/2) 


p{xi 


2/3) 


P(X2 


2/1) 


P(X2 


2/2) 


p{X2 


2/3) 


P{X3 


2/1) 


P{X3 


2/2) 


AX3 


2/3). 



(36) 



It is to be remembered that in our explicit example N = 
3. 



diag(a) = --diag((3 - a;^ - 2xf), (34a) 



diag(5) = --diag {{3~x^ - 2x) 



tanh 



diag(c) 



1 f^ + x 



-1-x 



1 



diag((3-y2_2y)2). 



diag(d) = -diag ((3 - y2 _ 2y) (1 + y)) 



(34b) 



(34c) 



(34d) 



diag(e) = diag {rm 



tanh 



-1 f^ + x 



0: 



tanh 



1 f^ + x" 



"tanh^i (i 



r)]' + [t-h-M^)]' 



.(34e) 



Replacing ipi with the polynomials in Eq. ( |12[ ), the col- 
location points of the problem then forms a mesh with 
the corresponding values p{xi^ yj) with i,j — 1, ..., A'^-l- 1. 
The collocation points are those illustrated in Fig. ([T]). 
However, instead of writing the polynomial as a matrix 



2. Boundary Condition Implementation 



Presently, Eq. (33) can be re-cast into matrix form as, 

Lp - Ap, (37) 

where the eigenvalues of the spectrum are obtained in A; 
the individual single particle energies. At this juncture, 
let us now suppose that we wish to solve the eigenvalue 
problem for the Isq state of the hydrogen atom. This 
state in the presence of a field would then be charac- 
terized using the notation i/^'^+^M'^^, where M = SiWi 
is the total z— component of angular momentum. The 
summation is over all the electrons in the atom; for the 
hydrogen atom there is only one electron. Similarly, the 
total z-parity of the state is the product of the z-parities 
of each electron, explicitly, tt^ — iTl^-^TTz^i- Again for 
hydrogen there is only one electron. The total z— com- 
ponent of the angular momentum, M, then forms a man- 
ifold within which different sub-spaces exist. The quan- 
tum number 1/ counts the excitation level within a given 
M— manifold and sub-space symmetry. The spin multi- 
plicity is given in the usual way as 2S + 1. Finally, the 
2;— parity of the state is indicated using tt^ = ±1, indi- 
cating positive or negative parity. Therefore, the state 
Iso of hydrogen in the presence of a strong field would 
be written as 1^0+. 

For this state of hydrogen the boundary conditions are 
as follows. Along both the x,y = +1 boundaries the 
wave function must vanish, therefore we have Dirichlet 
boundary conditions. Along the x,y = —1 boundaries 
however, we have Neumann conditions. These boundary 



conditions are to be kept in mind for the following discus- 
sion. A pictorial representation of the operator L acting 
upon the vector p is given in Fig. ([2]). 



N -I 



N -I 



N 



El 


Ci 


£■;, : Ci \ KV.I 




p(xi,yi) 

P{:':i,y2) 










p{xi,y3) 


£2 


C2 


El : Ci : £a'.2 




p{:':2,yi) 
p{:':2, J/2) 










P{^2,y3) 










p{=c3,yi) 
p(^3, y2) 

P(X3, i/s) 



N -1 



N -I 



N 



FIG. 2. Pictorial representation the operator L acting upon 
the vector p. The number of points in either direction, is 
A'' = 3 in the current example. 

Fig. (IT]) shows the different parts of the matrix operator 
L. Specifically, the different parts of the matrix that 
are relevant to the solution in the interior mesh points 
of the problem, are labelled. These are the submatrices 
Ei,E2, E3, E4, Ci, C2, C3, C4,En,i and En,2- The vector 
p can accordingly be split into three components, 



and 



P6„ 



Pbx 



'p{xi 


yiY 


p(xi 


y2) 


P{X2 


Vi) 


P[X2 


y2J 


Lp(a;2 


V3). 


p{xi 

P{X2 


,2/3) 
,2/3) 


P{X3, 
P{X3, 
P{X3, 


yi)" 

?/2) 
^3) 



(38a) 



(38b) 



(38c) 



In the above, Tpmt refers to the function value in the inte- 
rior points of the mesh given in Fig. (fTj) , while the func- 
tion values at the a; = — 1 boundary are given by Pf, 
and correspondingly, the values at the y = — 1 boundary 
are given by p^ . The eigenvalue problem in the interior 
mesh points is then given by. 



^Pint = '^Pint - C'Pb,, ~ EMVb 



(39) 



Compare Eq. (39 1 with Eq. (23). Meanwhile, the matri- 



ces E, C and E^ are given by(see Fig. ([2|)), 



E ■ 



El E3 

E2 E4 



(40a) 



C = 



Ci 
C2 



C3 

C4 



(7V-l)2x(7V-l) 



and 



E 



N 



N,l 

W,2 



(40b) 



(40c) 



{N-l)^xN 



The dimensions of the matrices are also quoted for con- 
venience. The objective here is to obtain a solution of 
the eigenvalue problem in Eq. ( 39 1 . This would be pos- 
sible if we are able to express the function values at the 
boundaries, i.e., p^ and Pf, in terms of the solution in 
the interior Pj^t. See discussion above regarding the one- 
dimensional problem; Eqs. (25-28). However, the values 
at the boundaries obey Neumann boundary conditions. 
A Neumann boundary condition imposed on the wave 
function ip along a given boundary is given by, 



n • VV' = g, 



(41) 



where, n is the unit vector normal to the boundary and 
g is the value to which, the directional derivative of the 
function along the direction of the normal vector is set. In 
our case, the boundaries in question are the lines x,y ^ 
— 1. The corresponding normal vectors are then trivial 
and we obtain the conditions. 



d_ 
dx 



ip{x,y) 



and 



dy 



ipix,y) 



(42) 



(43) 



y="i 



Once again, we can use differentiation matrices to specify 
these boundary conditions, as. 



b.pL=-i = 



and 



ByPl 



(44) 



(45) 



)Iy and 



The boundary matrices are defined by B^ — D^ 
By = Ix 'Si Dy. By virtue of the Kronecker product with 
the identity matrix, B^ has entries only along the diag- 
onals. The matrix By on the other hand, is a block di- 
agonal matrix with each block of dimension N x N. The 
products BxP and ByP are required by Eqs. (44) and 



(45) to vanish along the specified boundaries, therefore 



(Af-l)2x(Af-l)2 



we need only restrict our attention to those parts of the 
boundary matrices. Fig. ([3]) shows a pictorial represen- 
tation of the boundary matrix Bj. acting on the vector 
p. Since the derivative is required to vanish along the 
X — —I boundary, therefore we focus our attention on 
the part of the boundary matrix B^, that acts on the 
vector p along the boundary in question, i.e., the last N 
rows. Fig. (pi) shows the sub-matrices that are needed, 
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FIG. 3. Pictorial representation the operator B^ acting upon 
the vector p. The number of points in either direction, is 
N = S. 
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FIG. 4. Pictorial representation the operator By acting upon 
the vector p. The number of points in either direction, is 
TV = 3. 



these are labelled Bq 1, i?o,2, ^0,3, i?o,4 and Bjy. Thus, 
we can write Eq. (44) as. 



■SlPint + B2Pb^ + BNPb^ = 0. 



The matrices Bi and B2 in Eq. ( 46 ) are given by, 
and 



(46) 



(47a) 



B2 — [Bo,3 Bo^i\ 



0>4JjVx(Ar-i) 



(47b) 



The dimensions of the matrices are indicated for conve- 
nience. Eq. (46) has two unknowns, p^ and Pf, that 
need to be expressed in terms of Pj^j . Therefore we need 
another equation. This is provided by Eq. ( 45 ) . 



We 



have shown in Fig. (B a pictorial representation of the 
boundary matrix By acting upon the the vector p. The 
different sub-matrices are also shown in Fig. (B. We can 
then write down an equation for expressing the boundary 
condition in Eq. (45) as. 



GPint + Hpf,^ = 0. 
The matrices G and H are given by, 

G = 



and 



H = 



Gi G3 

G2 G4 


(Af-l)x(JV-l)2 


Hi Hj, 
H2 H4 


(Ar-l)x(Af-l) 



(48) 



(49a) 



(49b) 



It is to be kept in mind that since By is a block diagonal 
matrix, with each block of dimension iV x TV, only the 
block diagonal portion of G and the diagonal of H, are 



respectively non-zero. In the current explicit example, 
therefore we have G2 = G3 = [ ] and H2 = H3 = 0. 



With these definitions in place, we can employ Eqs. (46) 



and ( 48 ) to obtain expressions for the unknowns Pf, 
Pj, in terms of the values in the interior pjj^^ as, 



and 



Pb„ 



H-'Gp,, 



and 



Pb. 



B]^\Bi - B2H-'G)p,, 



(50) 



(51) 



Substituting Eqs. (50) and (51 1 in Eq. (39) we obtain an 



eigenvalue problem for the solution in the interior mesh 
points as. 



[E — E^Bj^ (i?i — B2H G) — 
Gi/^iG) Pi„t = Api„,. 



(52) 



Eq. ( 52 ) can then be solved as an eigenvalue problem 



using standard methods to obtain the eigenvalues A and 
the eigenvectors. See Section |IV] for details regarding the 
numerical method employed for solving the eigenvalue 
problem. For the moment, we turn our attention to ex- 
tending this methodology to the two-electron problem. 



B. The Two-Electron Problem 

The HF problem for the two electron atom can be writ- 
ten, using the matrix formalism detailed above in a com- 
pact form as, 



'Li + |diag[$ 



DM 



|diag[a£; 



-|diag[aB] L2 + |diag[$£,,2] 




^^1 



vV'2 



(53) 
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It is evident upon inspection that Eq. (53) is a coupled 



eigenvalue problem. The operators Li and L2 are the op- 
erators defined in Eq. ( 33 ) . The direct and exchange op- 
erators, diag[$£)^i] and diag[a£;] make the problem non- 
linear, as they depend upon the solutions ipi. However, 
the problem is linearized by estimating the direct and 
exchange interactions using wave functions from the pre- 
vious iteration. The exchange interaction still couples 
the two electrons and as such, we are still required to 
solve a coupled eigensystem. 

Once the hydrogenic problem has been solved to ob- 
tain initial estimates for the wave functions, these can 
then be employed for calculating the direct and exchange 
potentials. In order to continue with our explicit exam- 
ple, let us suppose that we wish to calculate the energy 
of the helium atom in the configuration l'^(— 1)"*", or in 
terms of field-free notation, lso2p_i. Thus the hydro- 
genic problem would first need to be solved for each of 
the two electrons in the configuration. Let us label the 
electrons' wave functions using different letters; p = Isq 
and q = 2p_i. 



The electron in the orbital 2p_i, in the presence of a 
magnetic field, has properties rather different from the 
Isq electron. The wave function of the former has posi- 
tive z— parity, but goes to zero along the y— axis. In other 
words, Q;, =0 and it has a Dirichlet condition imposed. 
There is however still a Neumann condition along the 
x-direction. This is to be remembered in the following 
discussion. Presently, we briefiy describe the method of 
solution of the elliptic partial differential equations for 
obtaining the direct and exchange interactions. 



1. The Direct and Exchange Interactions 

Let us assume that we have solved the hydrogenic prob- 
lem and already obtained initial estimates for the wave 
functions of each of the two electrons, viz., p and q for 
the states Isq and 2p_i, using the method described in 
Section III A Eqs. ^ and Q can be re-written after 



domain compactification as. 



i(3-.^-2.)^^ + i(3-a:^-2.) 



tanli 
1 



-1 f^ + x 



2(3-.^-2,)(l + ,)^^ 




$15,, (x,y) = -47r|V'j(a;,2;)| , (54) 



for the direct interaction and as 



i(3--^-2.f^ + ^(3~.^-2a:) 



tanh 



-1 A + a; 



2 (3 - y2 _ 2y) (1 + y) ^ - i^M ~ ^if {tanh-i (^-^^ 



for the exchange interaction. 



r 




"E {x, y) = -^ntp* (x, y) -ipi {x, y) , (55) 



More compactly, Eqs. (54) and (55) can be written using 
matrix form as, 



LdiY^D,i = -47r <, 2 



p^ for i = 1 
q^ for i == 2 



and 



ioxcha_B = -47r(pq). 



(56) 



(57) 



In the above, « = 1 or 2, labels the electrons. We empha- 



size at this stage that Eqs. (56 1 and (57) are written for 
the two-electron problem. For more than two electrons 
there would exist summations on right hand sides over 
the different states (see the direct and exchange terms in 



Eq. (jlj)). Moreover, p^, q^ and (pq) are the element-by- 
element products of the corresponding vectors that have 
be suitably normalised. Ldir and ioxch are the left hand 



side operators in Eqs. (54) and (55) respectively, the var- 



ious terms in which are defined similar to Eqs. (34a-e 



In the current explicit example, since p = Isq, the 
boundary conditions for the vector p^ are identical to 
that for p, i.e., Neumann conditions along x,y — —1. 
Likewise, since q = 2p_i, the boundary conditions for 
the vector q^ are identical to that for q, i.e., Neumann 
condition along y = —1 but Dirichlet condition along 
X = -1; qb^ = 0. 

The product of the two vectors (pq) on the other hand, 
vanishes along the line x = —1, since qi,^ =^ 0; therefore it 
has a Dirichlet condition given by (pq)6^ = 0. However, 
along the line y = — 1 a Neumann condition remains, 
since each individual vector has a Neumann condition. 

With these boundary conditions now identified, they 
can be imposed on the operators Ldir and iexch using 
the methods described in Section UlI A 21 The two linear 
systems of equations at the collocation points, given by 
Eqs. (56) and (57), are solved using standard methods. 



to obtain the direct and exchange potentials as, 



$D,i - {£;dir - CH-^G} ' qj 



2 
int ' 



(58a) 



^D,2 — {Edii - Ej^^dh -Bjv (-Bi - B2H G) 

CH-^Gy^-pt 



2 

int • 



and 



aE = {^exch - CH-^G] ' (pq) 



int ■ 



(58b) 



Be) 



In the above, the matrices -Ed 



E 



N,dir 



and -Eoxch are 

defined similarly to Eqs. (40a) and (c), this time how- 

and -Loxch respectively (see 
Section 



ever, using the operators Lji, 

Once the direct and exchange interac- 



IIIA2) 



tions have been determined, they can be substituted in 
Eq. ( 53 ) and the coupled eigensystem can be solved. 



We would like to caution at this stage, that the bound- 
ary conditions implemented in the explicit example are 
specifically for the configuration of the helium atom given 
by l'^(— 1)+ or lso2p_i. For other configurations, the 
boundary conditions imposed on p and q would be dif- 
ferent. In that case, Eqs. (58a-c) would change accord- 
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ingly. With this in mind, we now proceed to the next sec- 
tion which describes the setup of the coupled eigenvalue 
problem in Eq. ( 53 ) and the implementation of boundary 



conditions for its solution. 



2. The Coupled Eigenvalue Problem 

The direct and exchange potentials are found in 
Eqs. (|58afi 



-cj as vectors. These are converted to matri- 
ces with entries on the main diagonal before substituting 
into Eq. ( 53 ) . If we label the operator on the left hand 



side of Eq. (531 as M, then we can re-write Eq. (53 1 as. 



^Mn M12 
,M2i M22 




(59) 



We can depict pictorially, the action of Af on the vector 
, J or equivalently ) , as shown in Fig. ( 5 ) 



It can immediately be seen in Fig. ([5| that in the ma- 
trices Mil and M22, the off-diagonal sub-matrices are 
identical to those in Fig. (|3|, see Eqs. ( |40a| -c). Also, it 
can be seen that the matrices M12 and M21 are diagonal 
matrices that are identical. These represent the exchange 
interaction between the two electrons. Only the non-zero 
parts of these matrices that act on the interior parts of 
the vectors, Pi„t and qj„t are shown in Fig. ^. In this 
explicit example, the vector p represents the Isq state, 
while the vector q represents the 2p_i state. Accord- 



ingly, in the vector 



a Dirichlet boundary condition 



has been imposed explicitly along x = — 1 for the vector 
q by setting qj,^ — 0, as shown in Fig. ([s]). The coupled 
eigensystem can then be written as a system of coupled 
matrix equations for the interior points as. 



i?"Pi, 



T'\, 



APint - C'Pfc, - ENPb, (60a) 



and 



T 



12, 



E- 



'22, 



-^Qint - Cqi, 



(60b) 



The sub-m atrices G and En are defined as given in 



have slightly different entries on the diagonal and are 
thus defined as, 



E'' 



Ef E3 
E2 Ef 



1,2. 



(61) 



(Ar-l)2x(Ar-l)2 



I 

Similar to our discussion regarding the hydrogen atom in 
Section III A 2 we are required to express the vectors p;, , 
Pf, and q^ m terms of p^^^^ and qj^^ respectively, by im- 
plementing Neumann boundary conditions. This would 
enable us to then cast the coupled eigenvalue problem 
into its final form, prior to its numerical solution. 

Once again we construct the boundary matrices B^ 
and By as described earlier. This time however, we need 
to implement Neumann boundary conditions for both 
electrons. Thus, in addition to Eqs. (44) and (45), we 
have a Neumann boundary condition along y — —1 for 
the vector q, i.e., 

Byq\y^-i = 0. 



(62) 



Expressions for p;, , p;, and qj, , in terms of p-^^^ and qj^^^ 
respectively, are obtained using the method described 
earlier in Section IIII A2I These when substituted into 
Eqs. (|60a| and b) immediately yield. 



{£;" - CH-'G - EnB^\Bi - B2H-'G)} p,^^ 

^Pint 

(63a) 



+T'\n 



Eqs. (40b & c). However, The sub-matrix E^^ and E^'^ and 



Ti2p.„^ + {E'^ - CH-'G] qi„t = Aq, 



(63b) 



This can then immediately be recast as a standard cou- 
pled eigenvalue problem for the interior mesh points as. 
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FIG. 5. Pictorial representation the operator M acting upon 



the vector 
iV = 3. 



The number of points in either direction, is 



[S" - CH-'G - EnB],\Bi - B2H-^G)\ 



T 
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Eq. ( 64 1 is then solved using standard algorithms to ob- 
tain the eigenvectors and eigenvalues. Once the wave 
functions have been obtained, the direct and exchange 
interactions are re-evaluated. These updated potentials 
are then reused to carry out the HF iterations until con- 
vergence is achieved. We would like to caution the reader 
once more that the formulation of the HF problem shown 
in Eq. ( 64 1 is for the explicit example of the configura- 
tion of helium given by l'^(— 1)+ or lso2p_i. For other 
configurations, depending upon the boundary conditions, 
Eq. ( 64 1 will take a very different form. The left hand 



side operator shown in Eq. (64) is the pseudospectral 



representation of the HF operator for a particular con- 
figuration of the helium atom. 



[E"^^ -CH-^G]j \qi 



int > 



.qint, 



(64) 



It is also to be mentioned at this juncture, that we are 
not finding an approximate effective potential for the ex- 
change interaction as first suggested by Slater [7S] which 
greatly simplifies the HF equations. As a result Eq. (64) 



takes the form of a fully coupled eigenvalue problem, and 
not, we emphasize, the usual single-particle form, 

Fji/'j = ei/"*, (65) 

where F is the usual Fock operator given explicitly by 

F, = H, + Y.^Jj~K^). (66) 

i 

The direct (Jj) and exchange interactions {Kj) are usu- 
ally found using an approach similar to Slater's, employ- 
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ing wave functions from the previous iteration. This pro- 
cedure essentially uncouples the eigenvalue problem |75j . 
It can be seen upon examining Eq. (64) that we are not 
carrying out such an approximation. Therefore we solve 
the full HF problem as a coupled eigenvalue problem. 
In our computations for helium we found that the solu- 
tion of the uncoupled problem following Ref. [75] yielded 
slightly more bound results, by about 1 — 2%. That be- 
ing said, our goal here is to provide a fast and efficient 
means of solving the fully coupled problem, therefore the 
numerical procedure outlined below pertains as such, to 
the latter and the results presented in Section |V] are like- 
wise for a solution of the coupled problem. 

The method described above can easily be extended to 
tackle the case of the multi-electron atom, say lithium. 
It can be seen that the problem size in Eq. (64) will 



grow not only with the number of mesh points but also 
with the number of electrons. For a given number of 
mesh points N in each direction and a certain number of 
electrons rie, the size of the HF operator in Eq. (64 1 is 
[ne{N — 1)^) X {ne{N — 1)^] . If one were to express the 
exchange interactions completely using information from 
the previous iteration and thus uncouple the eigenvalue 
problem as in Eq. ( |65[ ), then within each HF iteration, 
one would be solving Ue eigenvalue problems each of size 
[{N — 1)^ X {N — 1)^] . Thus, since computer memory 
requirements are governed by this latter system size, the 
coupled problem in Eq. ([64]) is readily seen to be far more 
intensive than the uncoupled problem in Eq ( 65 ) . This 



important distinction is to be kept in mind when inter- 
preting the results obtained in this study. That being 
said, we give here presently, the block matrix form that 
the coupled HF problem will take for the lithium atom, 



/Mu Ti2 T3i\ /p\ /p\ 



Tl2 M22 T23 



Vt3i y23 mJ \rj 



(67) 



Vr/ 



The exchange operators are symmertric under permuta- 
tion, thus r*J = TJ\ i,j = 1, 2, 3. The matrices Mu, i = 
1,2,3, on the other hand, are formulated in a manner 
similar to that depicted pictorially in Fig. ([5]), modulo 
implementation of appropriate Dirichlet and Neumann 
boundary conditions, depending up the state in question. 
Thus for a large number of node points and electrons, 
it becomes imperative that efficient computational meth- 
ods be used in order to facilitate an economical solution 
of the fully coupled HF problem in Eq. (II]) , not only in 
terms of computer memory, but also quite importantly, 
in terms of computational time. The primary objective, 
as was stated earlier, is to delineate a method for the 
expedient solution of the single-configuration HF prob- 
lem, providing a software that could be easily integrated 
into atmosphere models for neutron stars. It is to be ac- 
knowledged at the very outset, that these results can be 
improved by using CI or MCHF methods. It is also to be 
acknowledged that since we are not using the usual ap- 



proximation for the exchange interaction, the computer 
memory requirements are quite a lot higher than the un- 
coupled problem. This fact should be kept in mind should 
the atomic structure software described herein be consid- 
ered for integration with atmosphere models. 

The following section details the efficient numerical 
methods employed for this solution and the convergence 
tests that are carried out. 



IV. NUMERICAL DETAILS 

Examination of Fig. ^ reveals that the matrix M, 
is largely a sparse matrix. Thus we can take advantage 
of this fact and employ sparse matrix methods for the 
method of solution of the coupled eigenvalue problem. 

An atomic structure software package based on the 
pseudospectral method outlined in Section |III| was de- 
veloped for this study. The code was written in the high 
level programming language MATLAB® making partic- 
ular use of its fast matrix manipulation algorithms. The 
eigenvalue problems described in Eqs. ([52]), ([64|) and ([67|) 



are solved by discretizing the equations and solving the 
resultant algebraic eigenvalue problem. The discretiza- 
tion is done using a standard Chebyshev-Lobatto spectral 
collocation method [72]. The coupled eigenvalue prob- 
lem in Eq. ([I]) is then solved using a sparse matrix gen- 
eralized eigensystem solver using the implicitly restarted 
Arnoldi method (IRAM) [76H78] . A numerical implemen- 
tation of this method is readily available in the software 
package ARPACK |,78j. The key advantage of employ- 
ing IRAM is that the memory storage requirements are 
far less than the original Arnoldi algorithm. Very briefly, 
it employs a typical Arnoldi factorization to generate an 
orthogonal basis for forming a Krylov subspace. The 
implicit restarting is closely related to the well under- 
stood implicitly shifted QR algorithm [TQ] [80] , where the 
idea is to restart the Arnoldi factorization with a vec- 
tor that is better pre-conditioned so that it can damp 
unwanted components from the eigenvector expansion. 
At this stage the reader is referred to Refs. [771 ESI [SI] 
for details on the implementation and other computa- 
tional aspects of ARPACK. The computational time is 
remarkably reduced due to implicit restarting, particu- 
larly for computing a few eigenvalues in a given part of 
the spectrum [7S] • This method was found to yield accu- 
rate results for the energy eigenvalues of the first few 
eigenstates of helium and lithium in intense magnetic 
fields. In our computations during a typical Arnoldi it- 
eration, it was seen that generating a Krylov subspace 
with about 50 basis vectors was sufficient for determin- 
ing around 15 eigenvalues in the vicinity of a given shift 
(cr), by employing the shift-invert algorithm [75] • Runs 
were carried out for different values of the magnetic field 
strength parameter j3z, in the range 5x 10~^ < Pz 1^ 10"^, 
for the cylindrical pseudospectral code. A typical toler- 
ance of around 10" ^'^ was employed for the internal errors 
of ARPACK. Additionally for testing convergence of the 
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pseudospectral method, we employed up to six different 
levels of m.esli refinem.ent, ranging from coarse to fine 
mesh i.e., N — 21,31,41,51,61 and 71 mesh points in 
each direction. The finest mesh calculations for lithium 
for example, took on the order of 1000 — 1400 seconds 
of computing time on an Intel® Xeon® E5620 2.4 GHz 
processor, for obtaining a fully converged HF solution. 

The code takes as its input, the number of electrons in 
the atom ?ie, the nuclear charge Z, and the magnetic 
field strength parameter (3z smd the configuration for 
which calculations are to be carried out, such as the state 
lso2p_i3d_2, etc. It then proceeds to compute system- 
atically the eigenvalues and eigenfunctions of the coupled 
system of equations in Eq. (II]) , according to the iterative 
procedure described below in brief. 

Eqs. [l] |3]and|4]are solved in a three step process us- 
ing the iterative self-consistent field (SCF) Hartree-Fock 
method [17]. First, an initial estimate is obtained for the 
eigenvectors by solving Eq. [l] without the contributions 
due the interaction between the electrons. The second 
step involves obtaining estimates for the potentials due 
to the direct and exchange interactions amongst the elec- 
trons, using the elliptic partial differential Eqs. [3] and |4] 
These are solved using the estimates for the wave func- 
tions obtained in the previous step. These potentials 
are then used to solve for better estimates of the eigen- 
functions along with the relevant eigenvalues in Eq. [T] 
The last two steps are iterated in the order described 
above to obtain progressively better estimates for the 
eigenvalues and eigenvectors with each HF iteration. It 
was observed during our runs that fast convergence was 
achieved; within about 6 — 12 HF iterations. A con- 
vergence criterion for the HF iterations was employed 
wherein the difference between the HF energies for two 
consecutive iterations was tested. Typically, a tolerance 
on the order of lO^^i^^^oo was employed for this purpose. 
Once the HF iterations attained convergence, the total 
energy of the Hartree-Fock state under consideration is 
reported according to Eq. [2] 

Additionally, the current work does not include rel- 
ativistic corrections to the energies. These corrections 
have been shown to be small in strong as well as intense 
magnetic fields for hydrogenic atoms IHMHl]- Their re- 
sults for the hydrogen atom reveal that for the states con- 
sidered in their studies, the fractional change in energy 
is on the order of 10~^ or so. This fractional change in 
energy was seen to be smaller than the numerical errors 
arising from convergence of the pseudospectral method, 
established using different levels of mesh refinements. 
Thus, while relativistic corrections are important, it was 
not possible to account for them accurately in the current 
study. 

Moreover, it is to be kept in mind that the current 
implementation of the pseudospectral method for the so- 
lution of the HF problem, is a single configuration cal- 
culation. Thus, it does not take into account effects 
such as electron correlation which can become impor- 
tant for many-electron systems. Post-HF methods, such 



as configuration-interaction and MCHF methods will no 
doubt yield more accuracy. However, as mentioned in 
Section |T] the aim of the current study is to facilitate a 
fast and computationally economical method for the 2D 
solution of the many-electron single-configuration prob- 
lem, without resorting to any basis expansion, or sepa- 
ration of variables or even the commonly employed adi- 
abatic approximation. The energies and wave functions 
obtained from such a calculation, could be directly em- 
ployed in a 2D-configuration-interaction calculation, or 
the method described in Section IIIII could be extended 
towards a full 2D-MCHF framework. Such an undertak- 
ing, while being perhaps the next step in the evolution 
of the atomic structure code described herein, was con- 
sidered to be outside the scope of the current paper. 

Presently, the results are presented alongside a discus- 
sion in the following section. 



V. RESULTS AND DISCUSSION 

The atomic structure software package developed in 
the current study extends our earlier computations for 
atoms in strong magnetic fields ^46, 48^ , towards the in- 
tense field regime, /3z :» 1. The states that were consid- 
ered in this study are labelled using both the field-free 
and strong-field notations for the sake of clarity; these 
can be found in Table U which lists these different states 
of helium and lithium. We have studied the six most 
tightly bound states of each atom in the intense magnetic 
field regime {Pz >• !)■ Within a given parity sub-space, 
typically there are crossovers that occur as the magnetic 
field is reduced and the reader is referred to Ref. |35] for 
an excellent discussion regarding ground state crossovers 
which are typically in the regime f5z < I. Within each 
given parity sub-space, we considered the most tightly 
bound state in the intense field regime. 



TABLE I. The different states of helium and lithium consid- 
ered in this study, listed using both intense-field and field-free 
notation. It is the the field-free configurations that are calcu- 
lated in the weak- and intermediate-field regimes. 





Intense-field 


Ficld-frcc 




1^(-1)+ 


lso2p_i 




l'(-l)- 


lso3d_i 


Helium 


l3(-2)+ 
l'(-2)- 


lso3d_2 

lSo4/-2 




1^(0)+ 


lso2so 




i'(o)- 


lso2po 




l*(-3)+ 


lso2p_i3ii~2 




l*(-3)- 


lso2p_i4/_2 


Lithium 


l4(-2)+ 
l4(-2)- 


lso2so3d_2 
lso2p_i3d-i 




l'(-l)+ 


lso2so2p_i 




1*(-1)- 


lso2po2p-i 
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A. The Helium Atom 

For the states of helium hsted in Table [ij eigenvalue 
computations were carried out and the total HF ener- 
gies were calculated for six different levels of mesh refine- 
ment. This enables us to extrapolate the convergence of 
the results to the limit of infinitely fine mesh, as shown in 
Fig. [6j This figure shows the convergence of eigenvalues 
obtained from a few sample calculations for both helium 
and lithium. With each level of mesh refinement, the sys- 
tem size grows and the full system size is reported on the 
a;— axis of Fig. |6] This is the dimension of the coupled 
eigenvalue problem that is solved. A variety of extrap- 
olation methods were tested and it was seen that using 
piece-wise continuous cubic Hermite interpolating poly- 
nomials [5S], yielded the most reasonably extrapolated 
results over the entire range of magnetic field strengths 
considered here. For extrapolation to infinitely fine mesh, 
we employed the average area per unit grid size as the 
abscissa. Thus it was possible to extrapolate to infinitely 
fine mesh or zero area per unit grid size. It is the ex- 
trapolated values that are reported in the tables for the 
binding energies. 



The calculated HF energies for positive parity states of 
helium are given in Table [ill while those for the negative 
parity states are shown in Table |III[ The absolute val- 
ues of the binding energies are given therein. Also shown 
therein, are the calculated data from other references for 
comparison. In this study we investigated the three most 
tightly bound states within each parity sub-space, in the 
limit of intense magnetic fields. The corresponding weak 
field orbitals of these states are those listed in Table [J 
The values given in parentheses are from eigenvalue com- 
putations using spherical coordinates. These calculations 
were carried out using our pseudospectral atomic struc- 
ture software developed in an earlier study [JS]. An im- 
proved and faster version of the code was employed, again 
using the same levels of mesh refinement for maintain- 
ing consistency. The cylindrical pseudospectral method 
begins to lose accuracy as the magnetic field decreases, 
in the weak field region, while in contrast, the spherical 
pseudospectral method loses accuracy in the intense field 
limit. Therefore using a combination of the two types of 
codes, we can explore the entire range in < /^^ < 1000. 



The positive z— parity sub-space of Helium 
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FIG. 6. Convergence of the binding energy with mesh re- 
finement is shown. The plot shows results from eigenvalue 
computations for two states of helium (top two panels) and 
lithium (bottom two panels) each, at different magnetic field 
strengths. The levels of mesh refinement employed correspond 
to A'^ = 21, 31, 41, 51, 61 and 71 points in each of the x— and 
y— directions for helium and up to A'^ = 61 for lithium. The 
lines drawn through the data are interpolating piece-wise cu- 
bic hermite polynomials that are also employed for extrapo- 
lation. 



It can be seen upon examining the data in Table |Tl] 
that the results obtained in the current study are in good 
agreement with values obtained elsewhere and in some in- 
stances represent slight to moderate improvements in the 
upper bound of the binding energy. For the state l'^(O)"'" 
which is comprised of the oribtals IsqIsq, we see that the 
average difference with respect to the most bound values 
given therein from other calculations, is about A « 1.25% 
with the maximum difference being at around j3z = 6.25 
at around 2%. The results remain accurate to large val- 
ues of magnetic field strength, i.e., f5z = 1000, where 
little data is available. The results obtained from the 
cylindrical pseudospectral method are about 0.8% more 
bound than the results of Ref. fB] at around Pz = 1000. 

Similarly for the state l'^(— 2)+, we notice that the av- 
erage difference from the most tightly bound values ob- 
tained elsewhere is about A sa 2.5%. In particular we 
see that both the cylindrical and spherical pseudospec- 
tral codes are slightly more accurate than our earlier 
study |46j in the intermediate and strong field regime; 
10~^ < /3z < 10, by a maximum of about 5%. The 
field-free spectroscopic orbital of the second electron in 
the state l'^(— 2)+, is the 3d_2 orbital, which has a much 
greater spatial extent than the Isq orbital, even in strong 
and intense magnetic fields. The same is true for the 2so 
electron in the configuration l''(0)+. In our earlier study 
[46] we had adopted a finite-element based approach for 
the solution of the single-configuration HF problem and 
as a result, had truncated the domain of the problem 
to about 20 Bohr radii in both the z— and p— direc- 
tions. However, in the current study we are employing 
a compactification which does not truncate the domain, 
but rather preserves the information in the entire two- 
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dimensional space. As a result, orbitals such 3d_2 that 
have a greater spatial extent, are determined with a little 
more accuracy in comparison to our previous study [15] , 
Concordantly, we see that the improvement in the strong 
field regime is on the order of 2 — 5% with respect to 
the results obtained in Ref. ^46j. We also see that in the 
limit of intense fields, i.e. Pz ^ 10, the accuracy is still 
preserved and we see an improvement in the estimates 
with respect to the Quantum Monte Carlo method due to 
Ref. [22] ■ While the method developed by the authors in 
Ref. [2HI11] may be more involved than the simple single- 
configuration HF treatment given here, nonetheless it is 
to be remembered that in the current study we are solv- 
ing a fully coupled eigenvalue problem for the electrons 
at each HF iteration and are not solving the problem 
by carrying out an approximation as suggested by Slater 
[75] (see discussion regarding Eqs. (65) and (66) above). 



This makes the current solution of the problem compu- 
tationally more intensive. Thus, overall in the current 
study, the computational overhead may be slightly more 
than that of the work in Refs. [22l [24]. With this cau- 
tionary note to the reader in place, we see that the pseu- 
dospectral code developed in this study produces slightly 
improved results. The IRAM method that is employed 



for calculating the eigenvalues (see Section IV) however. 



reduces the computational overhead quite significantly, 
in comparison to solving for the entire spectrum of the 
fully coupled problem. On average the computation of 
the converged values in Tables [ll] and |III| took between 
about 250 — 600 s of computing time for the finest mesh 
calculations, utilizing about 4 — 5GB of computer mem- 
ory 

The third state that was investigated was the ground 
state of helium in intense magnetic fields. This is the 
configuration 1^(— 1)"*" comprised of the spectroscopic or- 
bitals lso2p_i . Here we see that the values obtained from 
the pseudospectral method produce values that are in 
general a little more bound than the results obtained in 
our earlier study Ref. 46,. The maximum improvement 
in the range Q < Pz < 1.25 is about A w 0.7%. However 
we see that the configuration interaction calculations due 
to Schmelcher and co-workers, are in general more accu- 
rate than the single configuration results, although the 
results are still in reasonable agreement (A w 1.6% on 
average, over the range 2 < /32 < 12.5). This shows 
that taking into account the effects of electron correla- 
tion and interactions between different configurations is 
an important aspect necessary for the accurate solution 
of the multi-electron problem. The drawback however 
is that such calculations are computationally more in- 
volved than a simple single configuration calculation. We 
also see that for this state towards the higher end of the 
intense field regime, i.e., Pz ~ 1000, there is a minor 
loss in accuracy. This is attributed to the fact that the 
2p„i electron becomes exceedingly bound and concomi- 
tantly, greatly reduced in spatial extent and constrained 
closer to the nucleus while at the same time being elon- 
gated along the z— direction. Our finest mesh calculation 



yielded an absolute value for the HF binding energy of 
about 25.9923£'2_oo, however, the result obtained from 
extrapolation is what is quoted in Table [TTj Thus, we 
see that this result could be made more accurate by in- 
creasing the levels of mesh refinement so that better con- 
vergence could be achieved for the 2p_i orbital in the 
extreme end of the intense field regime. Our calculations 
were ultimately limited by computing resources, thus fur- 
ther mesh refinement in excess of A^ = 71 node points in 
each direction was not possible. 

Overall, for the l'^(— 1)+ state of helium the agreement 
with the most bound energy values from different stud- 
ies, is on average about A « 0.9%. Since our primary 
objective in this study is to delineate a fast and reason- 
ably accurate method for atomic structure calculations in 
intense magnetic fields, we considered this level of agree- 
ment to be sufficient. 

Next we discuss the three most tightly bound nega- 
tive parity states of the helium atom in intense magnetic 
fields. 



2. The negative z— parity sub-space of Helium 

The three most tightly bound negative parity states 
of helium in the intense field regime, were investigated in 
this stud y. T he results for the binding energies are shown 



in Table 



HI 



For the state l'^(O)^ we see that the results 
from the pseudospectral calculations, in both weak and 
strong fields {0 < (3z < 10), are in general more bound 
in comparison to the results obtained in our earlier study 
[46] . The average improvement is about A « 3.2%. As 
mentioned earlier, the current study employs a compact- 
ification of the entire domain without truncation and as 
a result the estimation of the orbitals that have a greater 
spatial extent is slightly more accurate in comparison to 
Ref. [in]- We also see a moderate improvement in the 
intense field regime where the average improvement is 
around A « 4% with respect to the most bound results 
quoted in the table. At the higher end of the intense 
field regime (/3z > 50), the results continue to remain 
accurate, as the pseudospectral method better estimates 
the spatial extent and shape of the 2po electron in the 
l'^(0)~. Here the improvements with respect to the data 
from Refs. [5] and [^ are about A « 2%. 

The calculated binding energies of the remaining two 
tightly bound states l'^(— 1)~ and 1^(— 2)~ show a sim- 
ilar trend. The average improvements with respect to 
the most bound results obtained elsewhere are around 
A sa 3.8% in both cases. It can be seen upon exam- 
ining the data for these states, that the state l'^(— 2)~ 
becomes more bound than the 1^( — 1)~ state of helium 
with increasing magnetic field strength in the intense field 
regime. However the binding energies are still quite close 
to each other. 

One of the aims of the current study is to provide a 
fast method for the calculation of the energy landscape of 
atoms in intense magnetic fields; therefore, we have addi- 
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tionally calculated fits to the data provided in Tables [TT] 
and mil The model fits are rational functions whose an- 
alytic form is given by, 



/(^) = 



En i 

i=0 ^i^ 






(68) 



where x = ln(l + /3) and m = n — 2. The fitting was 
carried out using a non-linear least squares Levenberg- 
Marquardt algorithm with line searches [86] . The coeffi- 
cients and the maximal fitting errors o ver the entire range 
/32 = to /3z = 10^ are given in Table IV As can be seen 



in Table IV the errors are small enough that these fitting 
functions could be employed directly in atmosphere mod- 
els of neutron stars rather than incorporating a code that 
calculates the binding energy. Thus, atmosphere models 
which are computationally intensive to begin with, need 
not be further complicated with the addition of an atomic 
structure calculation module, even though the software 
developed in this study is compact and computationally 
efficient and economical in terms of computational time 
as well. 



B. The Lithium Atom 

We investigated the six most tightly bound states of 
the lithium atom in intense magnetic fields. The bind- 
ing energies obtained from solving the coupled eigenvalue 
problem in Eq. (67) are shown in Tables [v| and VI These 
tables show the results for the positive and negative par- 
ity states, respectively. As in the case of the helium atom, 
the HF binding energies are results that were obtained 
after extrapolating to the limit of infinitely fine mesh. 



1. The positive z— parity sub-space of Lithium 

In contrast to the helium atom, lithium has been in- 
vestigated far less often in the literature, and data is 
scarce for the binding energies of the different states, 
particularly in the intense field regime. Table IV] shows 
the computed HF binding energies for the three most 
tightly bound states in the positive parity sub-space 
of the lithium atom. Once more, the values given in 
parentheses are those computed using our spherical pseu- 
dospectral code. Upon examining the data, we see that 
the state 1^(— S)"*" of lithium comprised of the orbitals 
lso2p-i3(i_2 is the best investigated state in the liter- 
ature thus far. The results obtained from the current 
study can be seen to be in reasonably good agreement 
with the estimates obtained elsewhere. The average im- 
provement in the estimate of the binding energy for the 
state l''(— 3)+, over the entire range < /^^ < 1000, is 
about A « 1%, with a maximum of around 3% in the 
strong field regime. The majority of the results for the 
binding energy of the 1*(— 3)"*" from other studies in Ta- 
blelvj are from Ref. [3B] . This study was also a single con- 
figuration calculation, albeit using a numerical approach 



employing finite-difference based mesh methods. We see 
that the current pseudospectral approach preserves ac- 
curacy as the magnetic field increases obtaining slightly 
better estimates than Ref. 36], even at the higher end 
of the intense field regime. Including effects of electron 
correlation and relativistic effects will of course improve 
the results obtained herein. 

The second most tightly bound positive parity state of 
lithium is 1"'(— 1) + . This state comprises of the orbitals 
lso2so2p-i. It can be seen that the estimates obtained 
in the current study are slight improvements with an av- 
erage of A « 0.8% over the range < /3z 5 56. While 
this state is not as tightly bound as the 1^(— 3)+ state 
in the intense field regime, examination of the data in 
weak fields reveals that with decreasing magnetic field 
strength, this state becomes the most tightly bound of 
the three states, in the vicinity of /3z ~ 0.3. The third 
state investigated is 1^(— 2)+, which represents the third 
most tightly bound state of the positive z— parity states 
in intense magnetic fields. For this state however, there 
is no data in the intense field regime. The only available 
data shown in Table [V] is in the weak and strong field 
regime. We see that in this case the average difference 
from the results obtained in Ref. [BS] is about A w 0.5%. 
And we see that the CI calculations produce a more ac- 
curate result at the higher end of the strong field regime, 
i.e., /3z « 0.5. To the best of our knowledge, since data is 
not available in the literature for this state in the range 
0.7 < Pz ^ 1000, no comparisons could be made and 
we state that the results obtained in the current study, 
represent the first such investigation. The same is true 
for the state l'^(-l)+ in the range 70 < Pz < 1000 



2. The negative z— parity sub-space of Lithium 



Table VI shows the three most tightly bound states of 
lithium that have negative z— parity. It is quite striking 
to see that the negative parity sub-space has not been 
tackled much in the literature. To the best of our knowl- 
edge, there is data available only for the !'*(— 1)^ state 
of lithium, comprised of the orbitals lso2po2p-i. We 
see that the pseudospectral approach produces estimates 
that are on average about A « 2.3% more bound than 
the results from other studies, over the range of magnetic 
field strengths < Pz ^ 56. However, once more there 
is no data available in the higher end of the intense field 
regime of 70 < /3z < 1000. 

The most tightly bound negative parity state of lithium 
is seen to be the state T'(— 3)^, which is comprised of 
the orbitals lso2p_i4/_2- It can be seen that this state 
becomes the most tightly bound negative parity state at 
around f3z ~ 20. Below this field strength, the l'*(— 1)~ is 
the most tightly bound of the three negative parity states 
of lithium shown in Table IVII To the best of our knowl- 
edge, this crossover has not been reported elsewhere in 
the literature. In addition, it can be seen by comparing 
the binding energies reported in Tables [V] and |VI[ that 
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the state 1^(— 3)~, is also the second most tightly bound 
state of lithium in intense magnetic fields. Moreover, the 
third state shown therein, the !'*(— 2)^ state, comprised 
of the orbitals lso2p_i3d_i, also has not been investi- 
gated in the literature. This latter state represents the 
third most tightly bound negative parity state of lithium 
in the intense field regime. 

Thus overall we see that for the six most tightly bound 
states of lithium in intense magnetic fields, two of the 
states have not been investigated earlier at all (!"'(— 3)~ 
and 1^(— 2)~) and a third state (l'*(— 2)+) has not been 
investigated in the intense field regime. Therefore the 
results presented here appear to be the first of such stud- 
ies. In addition we see that for the remaining three states 
that were investigated, the binding energies obtained in 
the current study show slight improvements relative to 
the estimates obtained elsewhere. 

Furthermore, for the sake of facilitating atmosphere 
and crustal models of neutron stars, we have also car- 
ried out rational function fits to the data. Once again 
these analytic forms can be implemented directly in such 
codes, thereby circumventing the need for atomic struc- 
ture calculations altogether. The rational functions have 
the same functional form as those described in Eq. ( 68 ) 
above. The coefficients of these rational functions are 
given in Table |VII| alongside estimates of the fitting er- 
rors. 

This concludes our discussion of the results obtained 
in this study. In the following section, we summarize the 
findings alongside a brief discussion of further avenues 
for investigation. 



VI. 



CONCLUSION 



We present below a brief summary of the paper and 
thereafter a short discussion of possible avenues for fur- 
ther work. 



A. Summary 

In the current study we have investigated low-Z atoms, 
helium and lithium in intense magnetic fields. A two- 
dimensional single-configuration Hartree-Fock method, 
as described in Ref. ^Bj, was adopted. A key feature of 
the method is that the potentials for the inter-electronic 
interactions are obtained as solutions to the elliptic par- 
tial differential equations as given in Eqs. (|3| and Q. 
The HF equations in Eq. (IT]) are solved using the self con- 
sistent field method. In the current study the HF equa- 
tions are solved as a fully coupled eigenvalue problem, 
without expressing the exchange interactions as an effec- 
tive single-particle potential |75j . Thus it was observed 
that the system size grew as ne{N — 1)^ x ne{N — 1)^, 
with Tie the number of electrons in the coupled problem 
and N the number of grid points in each direction. 



A pseudospectral approach was adopted for the nu- 
merical solution of the problem using cylindrical coordi- 
nates so as to facilitate calculations in the intense field 
regime. The resulting semi-infinite domain of the prob- 
lem was kept in its entirety and a suitable compactifica- 
tion was carried out. Thereafter, domain discretization 
was achieved using the commonly employed Chebyshev- 
Lobatto spectral collocation method. It was seen that the 
transformed equations had a coordinate singularity along 
the axis; a remnant of using cylindrical coordinates. The 
singularity was excised from the doma in, by a translation 
of about 6 — 10^^^; see Section III A 



Additionally, we outlined in great detail, the formula- 
tion of the pseudospectral problem with a description of 
the method for imposing boundary conditions. For this 
purpose we gave explicit examples of the pseudospectral 
implementation for the ground states of hydrogen and he- 
lium in intense magnetic fields. The method developed 
in the current study is a generalization of the implemen- 
tation due to the authors in Ref. [74 . 

The resulting discretized and coupled eigenvalue prob- 
lem problem was solved using standard sparse matrix 
methods. The software package ARPACK was employed 
for computing eigenvalues in the desired part of the 
spectrum, obtaining a handful of eigenvalues and eigen- 
vectors. A major advantage of the implicitly restarted 
Arnoldi method is the drastically reduced computational 
overhead and memory requirements, even for very large 
system sizes. As a result, we were able to obtain accu- 
rate eigenvalues and eigenvectors for helium and lithium 
in intense magnetic fields. 

The key enabling advantage of the psuedospectral ap- 
proach is the immensely reduced computational time re- 
quired for obtaining accurate results; on the order of 
about a thousand seconds. In addition, since we have 
adopted here an unrestricted two-dimensional approach 
to the problem |46j , it has the distinct advantage that we 
do not require a basis of functions to describe the wave 
functions. Thus the wave functions obtained in the cur- 
rent unrestricted 2D approach can be thought of in effect, 
as those arising from the superposition of a large num- 
ber of basis functions. Simultaneously, we also do not 
impose a separation of coordinate variables in the func- 
tional form of the individual electron wave functions and 
adhere instead to the natural symmetries of the problem, 
i.e., we maintain azimuthal symmetry. 

We presented data for the six most tightly bound states 
of the helium atom in intense magnetic fields, in Tables [TT| 
and |III[ These were seen to be consistent with findings 
elsewhere. Similarly we investigated the six most tightly 
bound states of the lithium atom as well. However, we 
found that the data in the literature to be rather scarce 
for lithium. As a result we could only compare our re- 
sults for four of the six states characterized in this study. 
We also presented, apparently for the first time, calcula- 
tions for the binding energies for the states !''(— 2)~ and 
1^(— 3)~ of lithium. We find that the the latter state is 
also the second most tightly bound state of lithium in 
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intense magnetic fields. 

The work described herein was motivated primarily by 
the need to have accurately determined upper bounds for 
the binding energies of atoms in intense magnetic fields 
employing a computationally straight-forward implemen- 
tation. As the atomic structure software developed here 
is small and computationally economical, it produces ac- 
curate results within a short amount of computing time. 
As a result, it can be incorporated directly into atmo- 
sphere and crustal models for neutron stars. However, 
while this may be desirable, it may present an additional 
layer of computational complexity. The user may wish 
to circumvent this by employing the rational function fits 
given in Tables |IV| and |VII| These analytic forms, model 
the data in the range Q < Pz < 1000 and thus may sim- 
plify atmosphere and crustal models considerably. 



B. Avenues for future work 

At this juncture it is prudent to be aware of the limi- 
tations of the atomic structure software. First, since the 
code developed herein is merely a prototype, we have 
not included effects such as electron correlation, rela- 
tivistic corrections, finite nuclear mass effects and elec- 
tron screening. Conceivably, these effects once included, 
will result in further improvements to the binding en- 
ergies. These additions would collectively represent one 
of the directions in which the current software could be 
improved. Secondly, the cylindrical pseudospectral code 
begins to lose accuracy in the intermediate range of mag- 
netic field strengths, Pz ~ 0.5. At this point it becomes 
necessary to carry out the computations using spherical 
coordinates. In the current study we are carrying out 
this switch manually, when we notice a drop in accuracy 
at low fields; this could be automated by keeping track 
of the change in the eigenvector estimates as the field 
is lowered gradually, in the weak and intermediate field 
regimes. 

In the current study the memory requirements for the 
code are quite lean, considering that we are solving a fully 
coupled eigensystem. As a result the current study can 
be extended to tackle the case of mid-Z atoms such as 
beryllium, boron, carbon, oxygen etc., for which binding 
energy data is even more scarce than that for lithium. 
This would be the obvious next step in the utilization 
of the current version of the software for determining 
atomic structure in intense magnetic fields. In this re- 
gard it may possible to enhance computational efficiency 
by rendering certain parts of the code to execute on a par- 
allel architecture. In particular, the determination of the 
eigenvalues within each HF iteration could be computed 
on a distributed system employing PARPACK; a par- 
allel implementation of the implicitly restarted Arnoldi 
method. This would have a considerable impact on the 
usage of memory as well, since array storage requirements 
would be distributed over numerous nodes. In addition, 
at present in the serial version of the code, the inter- 



actions between each pair of electrons are calculated one 
after the next. As a result the computational time for de- 
termining the interaction matrices depends cumulatively 
on the number of pair- wise interactions. Clearly this be- 
comes an issue with increasing number of electrons in 
the atom. This computation would directly benefit from 
calculating each pair-wise interaction on a different node 
in parallel, thereby reducing computational time and the 
rate would then be limited by the slowest pair-wise com- 
putation rather than the cumulative time for all the pair- 
wise computations. Concordantly, with parallelization it 
would be possible to increase mesh refinement particu- 
larly for atoms with ng > 3. In the current prototype of 
the software, we were restricted to A^ = 61 node points 
for lithium, due to an upper limit on the available com- 
puter memory. 

Moreover, the current version of the code is only accu- 
rate for the low-lying states. For example it was not pos- 
sible to obtain accurate estimates of excited states within 
a given symmetry sub-space, such as say the 2* (—3) + 
state of lithium etc. This is to be acknowledged as a 
drawback although it can be related directly to comput- 
ing resources. In this study we are generating a reason- 
ably sized Krylov subspace for determining a handful of 
eigenvectors of the fully coupled HF problem. Our pri- 
mary concern was with regard to memory usage, which 
can grow rapidly with system size. An obvious step would 
be to decouple the electrons in the HF iterations by calcu- 
lating an effective single particle exchange potential [7S] 
and solve the effective single particle HF equations as 
shown in Eq. ( 65 1 . This would reduce the computational 



memory requirements as the number of electrons in the 
problem are increased. By doing so, it would be possi- 
ble to determine further excited states in the spectrum, 
than what can be resolved in the coupled HF problem, 
for a given Krylov sub-space size and computer memory. 
Second, the mesh refinement can be increased further to 
resolve a greater part of the spectrum better. Once more 
this is limited quite obviously by computer memory. Al- 
though, it may be possible to have non-uniform grids in 
the x~ and y— directions. At the moment, the number 
of points in each direction is the same. This makes the 
formulation of the pseudospectral method somewhat sim- 
pler. However, with increasing magnetic field strength 
the electron orbitals become not only greatly reduced in 
spatial extent but also elongated along the direction of 
the magnetic field. Thus, by introducing greater number 
of points in the y— direction than in the x— direction it 
may be possible to re-formulate the problem in a more 
computationally efficient manner. 

Finally, the software described in the current study 
could be extended towards a full 2D configuration in- 
teraction or even a 2D-MCHF framework. These post- 
HF methods would lead to an immediate improvement 
in the estimates of the binding energies obtained herein. 
The inclusion of different configurations and incorporat- 
ing concomitantly, correlation effects between configura- 
tions, would produce better results. This would partic- 
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ularly be the case, in the weak and intermediate field 
regimes where correlation effects can contribute to the 
binding energy reasonably significantly. In either case, 
the method developed herein would be central to such 
enhancements and as such, the current study represents 
the very first implementation of a cylindrical pseudospec- 
tral method for atomic structure calculations in intense 
magnetic fields. 
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TABLE IV. Coefficients of the different rational functions for 
fitting tfie six states of helium discussed. The absolute maxi- 
mum fractional error of the eigenvalue relative to the fit from 
I3z ~ to /3z = W^ is reported in the variable e 



State 


Coefficients 


State 


Coefficients 










fflo = 


0.045473505 




fflo = 


-1.7217123 




ffli = 


1.3416631 




ffli = 


-17.814316 




02 = 


-1.9945187 




ffl2 = 


-19.974289 




aa = 


-9.5055929 


l'(-l) + 


ffl3 = 


0.70874301 


l«(-2)- 


04 = 


0.56822486 


lso2p_i 


04 = 


-0.94000973 


lso4f_2 


ffl5 = 


-0.48530257 




bo = 


1.6042624 




feo = 


-0.044558692 




bi = 


9.5985466 




6i = 


-1.1134528 




e = 


6 X 10^3 




62 = 

e = 


3.776243 
5 X 10"^ 




ao = 


-0.26130228 




fflo = 


-0.058636537 




ffli = 


-3.9622156 




ffli = 


-7.202045 




ffl2 = 


7.69114 




02 = 


2.1397541 




a-A = 


7.1549007 




as = 


75.211465 




04 = 


-24.249904 




04 = 


6.6163959 


l3(-2) + 


Clb = 


2.1137011 


1«(-1)- 


as = 


3.6429931 


lso3d_2 


fle = 


-0.99080421 


lso3d_i 


ae = 


-0.044520254 




bo = 


0.2528006 




60 = 


0.05689405 




bi = 


2.5631582 




61 = 


6.6097765 




&2 = 


-9.7214586 




62 = 


-17.233698 




&3 = 


9.8712035 




63 = 


-19.204331 




e = 


5 X 10-^ 




e = 


3 X 10-^ 




do = 


-0.21346351 




ao = 


0.69370335 




ffli = 


-7.0538773 




Oi = 


0.0079933198 




ffl2 = 


-10.924921 




02 = 


-5.9113942 


1^(0)+ 


ffl3 = 


0.35137528 


i«(o)- 


0-3 = 


0.31958673 


lso2so 


04 = 


-0.5304386 


lso2po 


04 = 


-0.39964524 




bo = 


0.1956651 




bo = 


-0.64093914 




bi = 


5.8292412 




61 = 


1.5039961 




e = 


4 X 10"^ 




£ = 


1 X 10"^ 
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TABLE V. Absolute value of the binding energies of the posi- 
tive parity states of lithium. Energies are in units of Rydberg 
energies in the Coulomb potential of nuclear charge Z = 3 for 
lithium. Accurate data from other work is also provided for 
comparison. (/3z = ■y/2Z'^). The values given in parenthe- 
ses are obtained from a faster version of our spherical atomic 
structure code developed earlier [48j . 





1^(- 


-2) + 


1* 


(-1)+ 




l*(-3)+ 


Pz 


Here 


Ref. [65J 


Here 


Elsewhere 


Here 


Elsewhere 





(1.1492) 


1.1491 


(1.1968) 


1.1926=' 


(1.1357) 


1.1427" [1.1299]'' 


0.00056 


(1.1541) 


1.1544 


(1.2024) 


1.1969=^ 


(1.1425) 


1.1487 


0.0028 


(1.1780) 


1.1720 


(1.2194) 


1.2121'= 


(1.1652) 


1.1663 


0.0056 


(1.1961) 


1.1901 


(1.2390) 


1.2334^ 


(1.1897) 


1.1869 


0.0111 


(1.2267) 


1.2203 


(1.2735) 


1.2674=^ 


(1.2324) 


1.2278 


0.0278 


(1.2964) 


1.2886 


(1.3530) 


1.3463" 


(1.3354) 


1.3294 


0.0556 


(1.3912) 


1.3930 


(1.4529) 


1.4432" 


(1.4699) 


1.4627 


0.5 


2.2993 




2.4094 




2.6782 




0.5556 


2.3759 


2.4145 


2.4885 


2.4246" 


2.7697 


2.6572" 


0.7 


2.5557 




2.6749 




2.9822 




1 


2.8657 




3.0006 




3.3547 




1.1111 


2.9664 




3.1063 


3.0432'' 


3.4721 


3.3695'' 


2 


3.6114 




3.7837 




4.2352 




2.7778 


4.0396 




4.2333 


4.1781'' 


4.7438 


4.6779'' 


5 


4.9451 




5.1836 




5.8237 




5.5556 


5.1282 




5.3756 


5.3304" 


6.0426 


6.0043'' 


7 


5.5536 




5.8219 




6.5525 




10 


6.2788 




6.5822 




7.4242 




11.1111 


6.5099 




6.8244 


6.7909" 


7.7027 


7.6856" 


20 


7.9535 




8.3368 




9.4492 




27.7778 


18.8840 




9.3109 


9.2936" 


10.5810 


10.5685'' 


50 


10.7991 




11.3187 




12.9213 




55.5556 


11.1787 




11.7205 


11.7000" 


13.3871 


13.3464'' 


70 


12.0662 




12.6336 




14.4618 




100 


13.5238 




14.1628 




16.2704 




111.1111 










16.8412 


16.7294'' 


200 


16.8116 




17.5969 




20.3551 




277.7778 


18.5917 




19.4541 




22.5789 


22.2774'' 


500 


22.1657 




23.1796 




27.0683 




555.5556 


22.8620 




23.9049 




27.9464 


27.4029'' 


700 


24.4514 




25.5599 




29.9557 




1000 


27.0771 




28.2923 




33.2904 





" Ref. [65] 
'' Ref. ^ 
" Ref. [59] 
•^ Ref. HZ] 
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TABLE VI. Absolute value of the binding energies of the pos- 
itive parity states of lithium. Energies are in units of Rydberg 
energies in the Coulomb potential of nuclear charge Z = 3 for 
lithium. Accurate data from other work is also provided for 
comparison. (/3z = ■y/2Z'^). The values given in parenthe- 
ses are obtained from a faster version of our spherical atomic 
structure code developed earlier [48] . 





l*(-2)- 


1^ 


(-1)- 


l"(-3)- 


Pz 


Here 


Here 


Elsewhere 


Here 





(1.1400) 


(1.1687) 


1.1652^^ 


(1.1306) 


0.00056 


(1.1457) 


(1.1735) 


1.1695'^ 


(1.1361) 


0.0028 


(1.1644) 


(1.1909) 


1.1865" 


(1.1556) 


0.0056 


(1.1827) 


(1.2110) 


1.2065" 


(1.1755) 


0.0111 


(1.2183) 


(1.2469) 


1.2417" 


(1.2104) 


0.0278 


(1.3073) 


(1.3354) 


1.3297" 


(1.2964) 


0.0556 


(1.4248) 


(1.4500) 


1.4463" 


(1.4113) 


0.5 


2.5144 


2.5901 




2.4866 


0.5556 


2.5985 


2.6712 


2.4842'' 


2.5718 


0.7 


2.7941 


2.8611 




2.7752 


1 


3.1302 


3.1899 




3.1167 


1.1111 


3.2387 


3.2961 


3.1035^ 


3.2265 


2 


3.9287 


3.9735 




3.9241 


2.7778 


4.3831 


4.4211 


4.2319^ 


4.3826 


5 


5.3384 


5.3651 




5.3445 


5.5556 


5.5309 


5.5577 


5.3767'^ 


5.5379 


7 


5.9775 


6.0017 




5.9885 


10 


6.7403 


6.7636 




6.7562 


11.1111 


6.9858 


7.0095 


6.8298'' 


7.0060 


20 


8.5011 


8.5128 




8.5180 


27.7778 


9.4756 


9.4858 


9.3242^ 


9.4947 


50 


11.4806 


11.4889 




11.5022 


55.5556 


11.8781 


11.8861 


11.7269^ 


11.8999 


70 


12.7929 


12.8008 




12.8156 


100 


14.3263 


14.3347 




14.3504 


200 


17.7625 


17.7758 




17.7921 


277.7778 


19.6188 


19.6369 




19.6531 


500 


23.3394 


23.3705 




23.3862 


555.5556 


24.0634 


24.0974 




24.1130 


700 


25.7156 


25.7561 




25.7715 


1000 


28.4424 


28.4943 




28.5116 



'^ Ref. [SH] 
^ Ref. [36] 
" Ref. m\ 
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TABLE VII. Coefficients of the different rational functions 
for fitting the six states of lithium discussed. The absolute 
maximum fractional error of the eigenvalue relative to the fit 
from Pz ~ to Pz = 10^ is reported in the variable e 



State 


Coefficients 


State 




Coefficients 




ao = 


-535.72875 






fflo = 


-9.3555515 




tti = 


-5798.3255 






ffli = 


-88.389164 




02 = 


-3697.7541 






ffl2 = 


-62.302697 


1^-3)+ 


aa = 


-75.172313 


l'(- 


"3)- 


0.3 = 


-0.17161294 


ls„2p_i3d_2 


a4 = 


-110.3002 


lso2p- 


-l4f-2 


04 = 


-2.0600694 




&o = 


468.35702 






bo = 


8.2232868 




&i = 


2070.8459 






bl = 


34.29776 




e — 


9 X 10"^ 






e = 


7 X 10-^ 












ao = 


1.0841938 




tto = 


-3.9323519 






ffli = 


-1.9251408 




tti = 


-42.850795 






ffl2 = 


-148.33918 




a2 ~ 


-39.165039 






ffls = 


-101.67185 


1^-1) + 


0.3 = 


0.42649409 


l\- 


-1)" 


04 = 


-1.2270097 


lso2so2p_i 


a4 = 


-1.4493673 


lso2po2p_i 


ffl5 = 


-3.1352421 




&o = 


3.2679311 






bo = 


-0.92575833 




&i = 


19.93254 






bl = 


7.0913811 




e = 


7 X 10"^ 






fe2 = 

e = 


60.317953 
3 X 10"^ 




«() = 


-3.5455821 






fflo = 


-17.90019903 




0.1 = 


-37.85524 






ffli = 


-25.29627281 




02 = 


-35.186156 






ffl2 = 


-661.06955428 


l4(-2) + 


as = 


0.52456655 


l'(- 


-2)- 


ffl3 = 


-1219.55537379 


lso2so3d-2 


a4 = 


-1.3412409 


lso2p- 


-i3d-i 


ffl4 = 


-139.85464269 




&() = 


3.0591966 






bo = 


429.60387534 




&i = 


18.498756 






bi = 


121.92364211 




e = 


9 X 10"^ 






e = 


6 X 10"'^ 



